1. Stats de base, modèles linéaires, GLM
Distribution statistique
Décrit les probabilités d’apparition des valeurs d’une distribution donnée.
Distribution normale
On parle aussi souvent d’une distribution Gaussienne.
\[ f(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2} \]
où \(\mu\) et \(\sigma\) sont des paramètres, respectivement la moyenne et l’écart-type.
Forme de la distribution normale
x<-seq(-5,15,by=0.01)
mu<-5
sigma<-2
y<-(1/(sigma*sqrt(2*pi))*exp((-1/2)*((x-mu)/sigma)^2))
plot(x,y,type="l")
Simuler une distribution normale
On peut utiliser la fonction rnorm pour simuler des
valeurs provenant d’une distribution normale.
set.seed(123)
n<-1000
x<-rnorm(n,mu,sigma)
hist(x,breaks=50)
Interprétation
v<-seq(-5,15,by=0.01)
hist(x,breaks=50,freq=FALSE,xlim=range(v))
curve(dnorm(x,mu,sigma),add=TRUE)
y<-dnorm(v,mu,sigma)
polygon(c(v[v>=6 & v<=8],8,6), c(y[v>=6 & v<=8],0,0),col=alpha("red",0.5),border=NA)
L’aire sous la courbe est de 1 et la proportion de la surface en
rouge représente la probabilité d’obtenir une valeur entre 6 et 8 avec
une moyenne de \(\mu = 5\) et un
écart-type de \(\sigma = 2\). On peut
estimer cette probabilité avec la fonction pnorm intégrée à
R. Ce type de fonction est intégré pour plusieurs distribution
statistique (voir ?Distributions).
Probabilité d’obtenir une valeur?
pnorm(5,mean=5,sd=2) # probabilité d'obtenir une valeur <= 5## [1] 0.5pnorm(-2,mean=5,sd=2) # probabilité d'obtenir une valeur <= -2## [1] 0.0002326291pnorm(4,mean=5,sd=2) # probabilité d'obtenir une valeur <= 4## [1] 0.3085375Probabilité d’obtenir une valeur entre 6 et 8
pnorm(8,5,2)-pnorm(6,5,2)## [1] 0.2417303Distribution d’une variable
\[ X \sim \mathcal{N}(\mu,\,\sigma^{2})\,\]
\(X\) suit une distribution normale avec une moyenne \(\mu\) et une variance \(\sigma^{2}\)
Erreur-type
\[ se = \frac{sd}{\sqrt{n}} \]
ou
\[ \sigma_{\bar{x}} = \frac{\sigma_{x}}{\sqrt{n}} \]
n<-100000
pop<-rnorm(n,100,20) # population fictive de n individus avec une moyenne de 100 et écart-type de 20
ech<-sample(pop,30) # échantillon aléatoire de 30 individus
mean(ech) # moyenne## [1] 100.5529sd(ech)/sqrt(30) # erreur-type## [1] 3.27161Simulation
Ici, on refait cet échantillonnage nreps fois et on
calcule l’écart-type des différentes moyennes obtenues.
nreps<-1000
ech1000<-sapply(1:nreps,function(i){
s<-sample(pop,30)
mean(s)
})
hist(ech1000)
Écart-type des valeurs simulées
sd(ech1000)## [1] 3.734432Comparons ceci à l’erreur-type calculé à partir de notre seul échantillon:
sd(ech)/sqrt(30) # erreur-type## [1] 3.27161Interprétation de l’erreur-type
L’écart-type de ces moyennes correspond à l’erreur-type estimée à partir d’un seul échantillon. En d’autres mots, l’erreur-type représente:
l’écart-type si des moyennes obtenues si on refaisait plusieurs fois le même échantillonnage
une mesure de précision dans l’estimation de la moyenne
Le concept d’erreur-type s’applique aussi à d’autres paramètres estimés et non seulement à la moyenne.
L’erreur-type d’une statistique ou d’un paramètre est une estimation de l’écart-type de sa distribution d’échantillonnage.
Intervalle de confiance
Reprenons notre population fictive pop qui a une moyenne
de 100 et un écart-type de 20. On se rappelle qu’un intervalle de
confiance pour une moyenne provenant d’un échantillon distribué
normalement (ou presque) se calcule de la façon suivante:
\[\bar{x} ± t_{dl,\alpha/2}\sigma_{\bar{x}}\]
On peut calculer cet intervalle de confiance avec R:
mean(ech)+qt(0.025,df=length(ech)-1,lower.tail=FALSE)*(sd(ech)/sqrt(length(ech)))*c(-1,1)## [1] 93.86168 107.24407Simulation
Maintenant, reprenons note population fictive et calculons à chaque fois cet intervalle de confiance. Ensuite, calculons quelle est la proportion des intervalles qui contiennent la véritable moyenne de 100.
nreps<-1000
ech1000<-lapply(1:nreps,function(i){
s<-sample(pop,30)
c(mean(s),mean(s)+qt(0.025,df=length(s)-1,lower.tail=FALSE)*(sd(s)/sqrt(length(s)))*c(-1,1))
})
ci<-as.data.frame(do.call("rbind",ech1000))
names(ci)<-c("mean","lowerCI","upperCI")
sum(100>=ci$lowerCI & 100<=ci$upperCI)/nreps## [1] 0.952Interprétation de l’intervalle de confiance
Un intervalle de confiance à 95% veut dire que 95% des intervalles de confiance qu’on obtiendrait si on refaisait le même échantillonnage contiendront la véritable moyenne de la population.
Plus généralement, on devrait parler d’un paramètre (moyenne, variance, pente, etc.) plutôt que d’une moyenne.
On entend souvent à tort que cela veut dire qu’il y a 95% de chances que le paramètre soit à l’intérieur des bornes. Or le paramètre est considéré comme étant fixe et la probabilité porte sur le fait que l’intervalle contienne la véritable valeur du paramètre ou pas.
Tout comme l’erreur-type, l’intervalle de confiance est une mesure de précision dans l’estimation d’un paramètre.
Valeur de p
Un exemple avec le test de \(t\) et une population fictive.
set.seed(1)
n<-10000
pop1<-rnorm(n,20,3) # population fictive de n individus avec une moyenne de 20
pop2<-rnorm(n,21,3) # population fictive de n individus avec une moyenne de 21
ech1<-sample(pop1,10) # échantillon de 10 individus de chaque population
ech2<-sample(pop2,10)
test<-t.test(ech1,ech2)
test##
## Welch Two Sample t-test
##
## data: ech1 and ech2
## t = -1.189, df = 14.686, p-value = 0.2533
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -4.030060 1.147318
## sample estimates:
## mean of x mean of y
## 19.22124 20.66261Interprétation de la Valeur de p
plot(0,0,xlim=c(-4,4),ylim=c(0,0.5),type="n")
curve(dt(x,test$parameter),add=TRUE)
v<-seq(-5,15,by=0.01)
y<-dt(v,test$parameter)
polygon(c(v[v<=test$statistic],test$statistic,min(v)), c(y[v<=test$statistic],0,0),col=alpha("red",0.5),border=NA)
polygon(c(v[v>=abs(test$statistic)],max(v),abs(test$statistic)), c(y[v>=abs(test$statistic)],0,0),col=alpha("red",0.5),border=NA)
abline(v=rep(test$statistic,2)*c(1,-1),lwd=2,lty=3,col="red")
Simulations
Reprenons le test de \(t\) précédent, mais cette fois, prenons nos échantillons dans la même population (pop1) à chaque fois.
nreps<-1000
pvalue<-sapply(1:nreps,function(i){
ech1<-sample(pop1,5)
ech2<-sample(pop1,5)
test<-t.test(ech1,ech2)
test$p.value#statistic
})
hist(pvalue,breaks=50)
abline(v=0.05,col="red",lwd=2,lty=3)
sum(pvalue<=0.05)/length(pvalue)## [1] 0.051Modèles linéaires simples
On a ici la formulation classique ou habituelle d’un modèle linéaire simple.
\[ y = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + ... + \beta_{n}x_{n} + \epsilon \]
\[ \epsilon \sim \mathcal{N}(0,\,\sigma^{2})\ \]
Autre formulation
\[ \mu = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + ... + \beta_{n}x_{n}\]
\[ y \sim \mathcal{N}(\mu,\,\sigma^{2})\ \]
où
\(y\) = observations
\(\mu\) = prédicteur linéaire
On a donc un modèle représentant le lien entre nos variables et on cherche à estimer les paramètres de ce modèles à l’aide de données.
Suppositions de bases
Quelles sont les suppositions de bases d’un modèle linéaire tel que celui présenté?
Normalité des résidus?
Homoscédasticité de la variance (la variance ne dépend pas de ou ne varie pas avec \(x_{n}\))
Linéarité des relations entre \(x_{n}\) et \(y\)?
Corrélation entre les \(x_{n}\)?
Indépendence entre les observations?
Simulations
Une des meilleures façons de s’assurer qu’on comprend comment fonctionne un modèle
Permet d’étudier le comportement d’un modèle avec des valeurs connues
Simulons un modèle
# nb d'observations
n<-75
# valeurs fictives des variables explicatives selon une distribution uniforme
x1<-runif(n,0,100)
x2<-runif(n,0,10)
x3<-runif(n,0,1)
# valeurs des paramètres beta
b0<-100
b1<-0.2
b2<-1
b3<--10
# prédicteur linéaire
mu<-b0+b1*x1+b2*x2+b3*x3
# observations à partir d'une distribution normale avec un écart-type de 5
y<-rnorm(n,mu,5)
# données
d<-data.frame(y,x1,x2,x3)Fonction lm
La fonction lm permet de faire des modèles linéaires
simples, ce qui inclue les régressions simples, les ANOVAs, ANCOVAs,
régression multiples, etc.
Remarquez l’utilisation de l’argument data où la
fonction ira chercher les différents éléments spécifiés dans la
formule.
m<-lm(y~x1+x2+x3,data=d)Coefficients
On peut appliquer plusieurs fonctions pour extraire les éléments d’un
modèle, en particulier un modèle effectué avec la fonction
lm.
coef(m)## (Intercept) x1 x2 x3
## 96.9666068 0.1882304 1.0717263 -4.0315717Sommaire
La fonction summary est la plus utile pour extraire
l’information d’un modèle.
summary(m)##
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.518 -2.927 0.206 2.739 12.382
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 96.96661 1.80048 53.856 < 2e-16 ***
## x1 0.18823 0.01909 9.858 6.23e-15 ***
## x2 1.07173 0.18827 5.693 2.62e-07 ***
## x3 -4.03157 1.92616 -2.093 0.0399 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.02 on 71 degrees of freedom
## Multiple R-squared: 0.6444, Adjusted R-squared: 0.6294
## F-statistic: 42.89 on 3 and 71 DF, p-value: 6.301e-16Vérification des suppositions de bases
On peut utiliser la fonction plot qui lorsqu’elle est
appliquée sur un objet de classe lm produira des graphiques
permettant de vérifier certaines suppositions de bases ou conditions
d’application.
par(mfrow=c(2,2))
plot(m)
Fonctions accessoires
On peut également créer facilement certains de ces graphiques avec d’autres fonctions.
par(mfrow=c(1,3))
plot(resid(m),fitted(m))
hist(resid(m))
qqnorm(resid(m))
qqline(resid(m))
Facteurs
# simulations d'un facteur
d$fac<-factor(sample(c("A","B","C"),nrow(d),replace=TRUE))
d$y<-d$y+2*as.integer(d$fac)-1
# modèle modifié avec le facteur
m<-lm(y~x1+x2+x3+fac,data=d)
summary(m)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + fac, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.3549 -3.0036 0.3166 2.6305 13.2288
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 97.2099 1.9870 48.924 < 2e-16 ***
## x1 0.1861 0.0190 9.794 1.09e-14 ***
## x2 1.0464 0.1854 5.645 3.39e-07 ***
## x3 -3.8057 1.8956 -2.008 0.0486 *
## facB 1.9366 1.4416 1.343 0.1836
## facC 6.5096 1.3744 4.736 1.12e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.932 on 69 degrees of freedom
## Multiple R-squared: 0.696, Adjusted R-squared: 0.674
## F-statistic: 31.6 on 5 and 69 DF, p-value: < 2.2e-16Comprendre son modèle avec des graphiques
La fonction predict sert à calculer la valeur prédite
par le modèle pour chaque observation. Par exemple, ceci permet de
mettre en relation les valeurs observées et les valeurs prédites. La
fonction predict dans ce cas calcule la valeur prédite par
le modèle pour chaque observation dans la base de donnée.
plot(d$y,predict(m))
Soumettre de nouvelle valeurs à predict
Plus souvent, on veut illustrer l’effet de nos variables et pour
cela, il faut soumettre des valeurs à la fonction predict.
Cette fonction est très utile lorsque l’on veut un maximum de contrôle
sur les valeurs à prédire.
x<-seq(min(d$x1),max(d$x1),length.out=50)
nd<-data.frame(x1=x,x2=mean(d$x2),x3=mean(d$x3),fac="B") # on fixe les autres variables à leur valeur moyenne
p<-predict(m,newdata=nd)
plot(d$x1,d$y) # observations
lines(x,p) # valeurs prédites
Package visreg
Le package visreg
permet d’illustrer facilement les prédictions d’un modèle (ou les effets
marginaux des différentes variables explicatives). Il utilise en
arrière-plan la fonction predict.
library(visreg)
par(mfrow=c(2,2))
visreg(m)
Package ggeffects
Le package ggeffects
permet également d’illustrer les prédictions d’un modèle, mais ce
package se base sur l’utilisation du package ggplot2 pour construire les
graphiques. Le principe demeure le même: en arrière-plan la fonction
predict est utilisée. On peut également combiner les
différents graphiques générés en utilisant le package patchwork et
sa fonction wrap_plots.
library(ggeffects)
library(patchwork)
g<-ggpredict(m)
p<-plot(g,add.data=TRUE)
wrap_plots(p)
Interactions
Dans vos mots, qu’est-ce qu’une interaction?
C’est lorsque l’effet d’une variable dépend des valeurs d’une autre variable (ou de plusieurs).
Un exemple
# modèle modifié avec le facteur
m<-lm(y~x1+x2+x3+x1*fac,data=d)
summary(m)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x1 * fac, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.0406 -1.9485 0.1316 2.1335 14.2374
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 98.50371 2.15978 45.608 < 2e-16 ***
## x1 0.14204 0.03151 4.508 2.7e-05 ***
## x2 1.15443 0.18388 6.278 2.9e-08 ***
## x3 -2.44076 1.86546 -1.308 0.19521
## facB 0.66407 2.63788 0.252 0.80201
## facC -1.01455 3.01397 -0.337 0.73746
## x1:facB 0.01671 0.04332 0.386 0.70090
## x1:facC 0.13612 0.04871 2.794 0.00678 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.692 on 67 degrees of freedom
## Multiple R-squared: 0.7329, Adjusted R-squared: 0.705
## F-statistic: 26.26 on 7 and 67 DF, p-value: < 2.2e-16Illustration graphique
visreg(m,"x1",by="fac",overlay=TRUE)
On a une interaction significative, mais est-elle importante d’un point de vue biologique? Peut-être pas… Toujours distinguer entre la signification statistique et la signification biologique.
Une valeur de p significative c’est bien beau, mais ce qui importe vraiment c’est la taille de l’effet.
Interactions avec ggeffects
g<-ggpredict(m,terms=c("x1","fac"),message=FALSE)
plot(g,add.data=TRUE,jitter=FALSE)
Interprétion
Beaucoup plus facile avec un graphique qu’en regardant une table de coefficients.
En présence d’une interaction (significative), les effets simples sont plus difficilement interprétables et ne peuvent être interprétés sans faire référence à l’interaction (à moins d’ajustements particuliers voir Schielzeth 2010).
Si pour deux variables impliquées dans une interaction les effets simples ne sont pas significatifs, mais que leur interaction l’est, il faut considérer que ces deux variables ont un effet même si les coefficients associés aux effets simples ne sont pas significatifs.
Types d’erreurs
- Erreur de type I: rejeter l’hypothèse nulle, alors qu’elle est vraie (faux positif).
- Erreur de type II: ne pas rejeter l’hypothèse nulle, alors qu’elle est fausse (faux négatif).
Simulations
Simulons un modèle dans lequel les variables n’ont aucun effet.
Que se passera-t-il si on répète plusieurs fois ce scénario?
À quoi ressembleront nos valeurs de p?
Créons d’abord une fonction sim_model pour générer un
modèle linéaire fictif.
sim_model<-function(b0,b1,b2,b3){
n<-200
x1<-runif(n,0,100)
x2<-runif(n,0,10)
x3<-runif(n,0,1)
mu<-b0+b1*x1+b2*x2+b3*x3
y<-rnorm(n,mu,5)
d<-data.frame(y,x1,x2,x3)
m<-lm(y~x1+x2+x3,data=d)
m
}Simulons
Ensuite, simulons nsims modèles fictifs et emmagasinons
ces modèles dans une liste nommée list_models.
nsims<-500
list_models<-lapply(1:nsims,function(i){
sim_model(b0=100,b1=0,b2=0,b3=0)
})Si on illustre les valeurs de p obtenues pour chacun des coefficients associés aux variables \(x_{n}\) à l’aide d’un histogramme, à quoi ressemblera cet histogramme?
Valeurs de p obtenus
p<-lapply(list_models,function(i){
summary(i)$coef[2:4,4]
})
hist(unlist(p),breaks=seq(0,1,by=0.05),xlab="Valeurs de p")
Conclusions
Dans une certaine proportion des cas ( ~ 5% ), on conclut à un effet, alors qu’il n’y en a pas pas! C’est le fameux seuil \(\alpha\) de 0.05.
Sous l’hypothèse nulle, les valeurs de p sont distribuées uniformément entre 0 et 1.
Modèles linéaires standards
\[y \sim \mathcal{N}(\mu, \sigma^2)\]
\[\mu=\beta_0 + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_k X_k\]
Formulation matricielle
\[y \sim \mathcal{N}(\beta\mathbf{X}, \mathbf{I}\sigma^2)\]
\(\mathbf{X}\) est une matrice avec les valeurs des variables explicatives.
\(\mathbf{I}\) est une matrice identité qui a pour effet d’assigner une variance égale \(\sigma^2\) à toutes les observations.
Modèles linéaires généraux
\[y \sim \mathcal{N}(\beta\mathbf{X}, \mathbf{\Sigma})\]
\(\mathbf{\Sigma}\) est une matrice de variance-covariance qui permet beaucoup plus de flexibilité en terme de variance et de dépendance entre les observations. En général, des paramètres seront associés à la matrice \(\mathbf{\Sigma}\) pour permettre différentes situations. Dans le cas d’un modèle linéaire standard, on a \(\mathbf{\Sigma} = \mathbf{I}\sigma^2\).
Ex.: variance différente pour les niveaux d’un facteur, variance augmentant avec certaines valeurs de \(x\), autocorrélation temporelle ou spatiale, etc.
Ces différentes situations peuvent être obtenues entre autres avec
les packages nlme
(fonctions gls ou lme) ou glmmTMB.
Modèles linéaires généralisés (GLM)
\[y \sim \mathcal{ExpFam}(\theta, ...)\]
\[\mathcal{g}(\mu)=\beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k\]
Structure d’un GLM
\[y \sim \mathcal{ExpFam}(\theta, ...)\]
\[ \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k\]
\[ g(E(y)) = g(\mu) = \eta \]
Composante systématique
\[ \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k\]
Le prédicteur linéaire \(\eta\) permet de relier les variables explicatives à la variable réponse.
Composante aléatoire
\[ y \sim \mathcal{ExpFam}(\theta) \]
Les observations \(y\) suivent une distribution particulière appartenant à la famille des distributions exponentielles avec des paramètres \(\theta\). La distributions normale, de Poisson, Binomiale et Gamma font parties de la famille des distributions exponentielles.
Fonction de lien
\[ g(E(y)) = g(\mu) = \eta \]
La fonction de lien (link function) est une transformation qui permet de relier la réponse attendue \(E(y)\) (i.e. la moyenne) au prédicteur linéaire.
C’est ce qui permet de relier les variables explicatives à la variable réponse. Pour chaque distribution qui pourrait être utilisée dans le cadre d’un GLM, une fonction de lien dite canonique est utilisée par défaut, mais il existe plusieurs fonctions de liens possibles.
GLM Poisson
\[ y \sim \mathcal{Pois}(\lambda) \]
\[ log(\lambda) = \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]
On a donc des observations \(y\) qui suivent une distribution de Poisson avec un paramètre \(\lambda\). Une fonction de lien \(log\) permet de relier linéairement les variables explicatives à la valeur attendue \(E(y) = \lambda\) de la variable réponse.
Distribution de Poisson
Cette distribution est souvent utilisée pour décrire le nombre de fois qu’un événement survient dans un intervalle de temps ou d’espace.
\[f(k,\lambda) = Pr(y = k) = \frac{\lambda^{k}e^{-\lambda}}{k!}\]
Cette fonction permet de calculer la probabilité d’observer un événement \(k\) fois avec une distribution de Poisson ayant un paramètre \(\lambda\). La moyenne de cette distribution est égale à \(\lambda\).
Cette distribution a la propriété que \(\lambda = E(y) = Var(y)\). Autrement dit, la moyenne et la variance sont égales et correspondent au paramètre \(\lambda\) de la distribution.
D’où provient-elle?
Simulons une distribution de Poisson en distribuant aléatoirement des observations dans un carré 10 x 10. Ensuite, divisons ce carré en 100 cellules.
library(sf)
n<-200
x<-runif(n,0,10)
y<-runif(n,0,10)
sfc <- st_sfc(st_polygon(list(rbind(c(0,0), c(10,0), c(10,10), c(0,0)))))
grid <- st_make_grid(sfc, cellsize = 1)
plot(grid,axes=TRUE)
points(x,y)
Nb d’observations par cellule
Ensuite, calculons le nombre d’observations par cellule et illustrons la distribution du nombre d’observations pour chaque cellule. En moyenne, on devrait retrouver 2 observations par cellule (200 obs. / 100 cell. = 2).
o<-st_intersects(grid,st_as_sf(data.frame(x,y),coords=c("x","y")))
co<-sapply(o,length)
hist(co,breaks=seq(-0.5,max(co)+0.5,by=1),freq=FALSE,main="",xlab="Nb d'observations")
lines(0:max(co),dpois(0:max(co),lambda=2),type="b",cex=2)
legend("topright",lty=c(NA,1),legend=c("Observée","Théorique"),bty="n",cex=1,pch=c(22,21),pt.bg=c("gray","white"))
Ce serait la même chose si on simulait des observations aléatoires dans le temps et si on comptait le nb d’observations par intervalles de temps.
Simulons un modèle
Simulons un modèle à partir d’une distribution de Poisson.
n<-100
x<-runif(n,0,10)
lambda<-exp(-2+0.4*x)
y<-rpois(n,lambda)
d<-data.frame(y,x)
m<-glm(y~x,family="poisson",data=d)
summary(m)##
## Call:
## glm(formula = y ~ x, family = "poisson", data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.46287 0.30996 -7.946 1.93e-15 ***
## x 0.45294 0.03784 11.971 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 289.290 on 99 degrees of freedom
## Residual deviance: 85.894 on 98 degrees of freedom
## AIC: 245.83
##
## Number of Fisher Scoring iterations: 5Coefficients
Les coefficients estimés sont très similaires à ce qui a été utilisé pour générer le modèle (-2 et 0.4).
coef(m)## (Intercept) x
## -2.4628709 0.4529379Prédictions
Illustrons l’effet de la variable explicative.
visreg(m)
Sous quelle échelle?
visreg(m,scale="response")
Échelle du prédicteur linéaire ou de la réponse
par(mfrow=c(1,2))
visreg(m)
visreg(m,scale="response")
Remarquez que dans le graphique de gauche, les valeurs en y sont négatives ou positives. En fait, la fonction log permet de transformer les comptes allant de 0 à l’\(+\infty\) en valeurs allant de \(-\infty\) à \(+\infty\). Sous l’échelle log, les valeurs entre 0 et 1 sont négatives et les valeurs supérieures à 1 sont positives.
\[ log(1) = 0 \] \[ log(0.1) = -2.302586 \] \[ log(5) = 1.609438 \]
ggeffects
g<-ggpredict(m,terms=c("x [n=50]"))
plot(g,add.data=TRUE,jitter=FALSE)
predict
v<-seq(min(x),max(x),length.out=50)
plot(y~x,data=d)
p<-predict(m,data.frame(x=v))
lines(v,p)
Argument type
Par défaut, l’argument type = "link" et il faut
spécifier type = "response" pour obtenir les prédictions
sous l’échelle de la réponse.
v<-seq(min(x),max(x),length.out=50)
plot(y~x,data=d)
p<-predict(m,data.frame(x=v),type="response")
lines(v,p)
Suppositions de base
On pourrait être tenté d’utiliser la fonction plot pour
vérifier les conditions d’applications.
par(mfrow=c(2,2),mar=c(4,4,2,2))
plot(m)
Quelles sont les suppositions de bases?
\[ y \sim \mathcal{Pois}(\lambda) \]
\[ log(\lambda) = \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]
Certaines suppositions sont les mêmes que pour un modèle linéaire standard avec une distribution, mais d’autres diffèrent.
Vérification avec DHARMa
Une façon de vérifier les supposition de bases des modèles est d’utliser le modèle pour simuler des observations et ensuite comparer ces observations simulées aux observations réelles. Si les observations réelles ne ressemblent pas aux observations simulées à partir du modèle, c’est que le modèle ne représente pas bien ce qui se passe avec les données réelles et donc que qqch devrait être modifié dans le modèle.
Le package DHARMa se base sur ces simulations pour comparer les observations attendues par le modèle aux observations réelles. Plus spécifiquement, la méthode utilise les résidus quantiles pour comparer chaque observation aux observations attendues pour cette observations. Les résidus quantiles sont rapportés entre 0 et 1 et si le modèle est approprié pour les données, on devrait s’attendre à ce que les résidus soient distributés uniformément entre 0 et 1 pour toutes les observations.
library(DHARMa)
s<-simulateResiduals(m)
plot(s)
Si le modèle est appropré pour les données, on devrait voir ici une distribution uniforme des points et il ne devrait pas y avoir de patrons particuliers dans la distribution des points dans le graphique de droite. En particulier, les lignes de prédictions illustrées devraient être alignées sur les trois lignes représentant les quantiles 25, 50 et 75%. Dans le graphique quantile-quantile de gauche, les points devraient être alignés sur la droite.
Autres vérifications
Les problèmes de sous- ou sur- dispersion et de surabondance de 0 (zéro-inflation) sont très fréquents avec les GLM. Ces deux fonctions de DHARMa permettent d’évaluer si les observations réelles diffèrent des observations simulées en terme de présence de 0 ou de disperison. En rouge, ce sont les observations réelles et en gris les observations simulées. Si les observations réelles sont bien en dehors des observations simulées, C’est qu’il y a un problème avec le modèle et il faut un ajustment.
par(mfrow=c(1,2))
testDispersion(s)
testZeroInflation(s)
GLM binomial
\[ y \sim \mathcal{Binom}(n,p) \]
\[ logit(p) = log\left(\frac{p}{1-p}\right) = \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]
Régression logistique
\[ y \sim \mathcal{Binom}(1,p) = \mathcal{Bernoulli}(p)\]
La distribution de Bernoulli est un cas particulier de la distribution binomiale avec un tirage (e.g. un individu est vivant ou mort, un individu est mâle ou femelle, etc.)
Distribution binomiale
\[f(k,n,p) = Pr(y = k) = \binom{n}{k}p^k(1-p)^{n-k}\]
où
\[\binom{n}{k}=\frac{n!}{k!(n-k)!}\]
Ce qui représente le nombre de façons de prendre \(k\) éléments parmi \(n\) éléments.
Cette fonction permet de calculer la probabilité d’observer \(k\) succès si \(n\) tirages sont faits et que la probabilité de succès est de \(p\)
Cette distribution a la propriété que \(E(y) = np\) et \(Var(y) = np(1-p)\) , autrement dit, la moyenne et la variance sont entièrement détemrinées par la propbabilité de succès et par le nombre de tirages. On a donc une autre distribution assez restrictive.
Exemple
Supposons qu’on tire pile (0) ou face (1) 10 fois et qu’on calcule le nombre de face et qu’on répète cette exercice 500 fois. Quelle sera la distribution du nombre de faces obtenus?
n<-500
co<-sapply(1:n,function(i){
s<-sample(0:1,10,prob=c(0.5,0.5),replace=TRUE)
sum(s) # on additione les faces (1)
})
hist(co,breaks=seq(-0.5,max(co)+0.5,by=1),freq=FALSE,main="",xlab="Nb de faces")
lines(0:max(co),dbinom(0:max(co),prob=0.5,size=10),type="b",cex=2)
legend("topright",lty=c(NA,1),legend=c("Observée","Théorique"),bty="n",cex=1,pch=c(22,21),pt.bg=c("gray","white"))
Simulons un modèle
logit<-function(p){
log(p/(1-p))
}
inv.logit<-function(x){
exp(x)/(1+exp(x))
}
n<-200
x<-runif(n,0,10)
z<--2+0.5*x
y<-rbinom(n,size=1,prob=inv.logit(z))
d<-data.frame(y,x)
m<-glm(y~x,family="binomial",data=d)
summary(m)##
## Call:
## glm(formula = y ~ x, family = "binomial", data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.68335 0.43063 -6.231 4.63e-10 ***
## x 0.63864 0.08853 7.214 5.45e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 272.74 on 199 degrees of freedom
## Residual deviance: 185.42 on 198 degrees of freedom
## AIC: 189.42
##
## Number of Fisher Scoring iterations: 5Visualisation
Avec visreg
par(mfrow=c(1,2))
visreg(m)
visreg(m,scale="response")
Remarquez que dans le graphique de gauche, les valeurs en y sont négatives ou positives. En fait, la fonction logit permet de transformer les probabilités allant de 0 à 1 en valeurs allant de \(-\infty\) à \(+\infty\). Sous l’échelle logit, une valeur de 0 corrrespond à une probabilité de 0.5.
\[ log\left(\frac{p}{1-p}\right) = log\left(\frac{0.5}{1-0.5}\right) = log(1) = 0 \]
Visualisation
Avec ggeffects
g<-ggpredict(m,terms=c("x [n=50]"))
plot(g,add.data=TRUE)
Suppositions de base
par(mfrow=c(2,2))
s<-simulateResiduals(m)
testUniformity(s)
testQuantiles(s)
testDispersion(s)
testZeroInflation(s)
Problèmes
Les problèmes les plus fréquents sont la surdispersion et la surabondance de zéros (zero-inflation)
Surdispersion
C’est lorsque les données sont plus variables que ce qui est attendu par le modèle. Les distribution de Poisson et binomiale sont très restrictives par rapport à leur variance et les problèmes de surdispersion sont très fréquents.
En général, une surdisperion ignorée fera en sorte que l’incertitude sera sous-estimée (trop de confiance dans nos résultats).
Plusieurs causes possibles (variables explicatives manquantes, phénomènes d’aggrégation, hétérogénéité, etc.)
Solutions:
Poisson -> Négative binomiale (glm.nb,
glmmTMB, brms, etc.)
Binomiale -> Béta binomiale (glmmTMB,
brms, etc.)
Obervation-Level Random Effects (OLRE) (voir Harrison 2015 et Harrison 2014).
Exemple
Voici un exemple de GLM de Poisson avec de la surdispersion. Pour créer cette surdispersion, j’utilise une distribution négative binomiale pour générer les observations au lieu d’une distribution de Poisson.
n<-100
x<-runif(n,0,10)
lambda<-exp(-2+0.4*x)
y<-rnbinom(n,mu=lambda,size=0.5) #
d<-data.frame(y,x)
m<-glm(y~x,family="poisson",data=d)
g<-ggpredict(m,terms="x")
plot(g,add.data=TRUE)
Vérification avec DHARMa
Les résidus ne sont clairement pas distributés uniformément sur le second graphique et la surdispersion est évidente sur le 3e graphique.
s<-simulateResiduals(m)
par(mfrow=c(1,4))
testUniformity(s)
testQuantiles(s)
testDispersion(s)
testZeroInflation(s)
Correction avec glmmTMB
library(glmmTMB)
m<-glmmTMB(y~x,family="nbinom1",data=d)
s<-simulateResiduals(m)
par(mfrow=c(1,4))
testUniformity(s)
testQuantiles(s)
testDispersion(s)
testZeroInflation(s)
Sousdisperion
C’est lorsque les données sont moins variables que ce qui est attendu par le modèle. Ceci est beaucoup plus rare et probablement moins dommageable (e.g. taille de couvée/portée chez les oiseaux/mammifères voir Brooks et al. 2019).
En général, une sous-dispersion ignorée fera en sorte que l’incertitude sera surestimée (trop conservateurs dans nos résultats) ce qui est un moindre mal.
Solutions:
Poisson -> glmmTMB, family = "compois"
Zéro-inflation
Surabondance de 0 dans les observations par rapport à ce qui est attendu par le modèle. Il est possible aussi d’avoir une sous-abondance de 0 et/ou encore une surabondance de 1, mais ce dernier cas est un peu plus rare.
Solutions:
modèle zero-inflated
(glmmTMB,pscl,brms, etc.)
À noter qu’une surdipersion peut être due à une surabondance de 0 et vice-versa. Ce sont deux problèmes qui peuvent interagir.
Extensions
Gamma \(]0,\infty\)
Variable réponse continue, strictement positive et avec une asymétrie à droite.
Souvent une alternative à une transformation \(log\).
Ex.: concentration d’un composé, biomasse, etc.
Tweedie \([0,\infty\)
Variable réponse continue positive et incluant des 0.
Ex.: biomasse capturée, précipitations, etc.
Voir Lecomte et al. 2013
Beta \(]0,1[\) ou \([0,1]\)
Variable continue entre 0 et 1 exclusivement.
Pour proportion ou % ne provenant pas d’un ratio entre entiers.
Si 0 et 1 possibles, certaines transformations peuvent être utilisées.
Ex.: proportion d’une surface affectée, proportion de temps passé en alimentation, etc.
Voir Douma & Weedon 2019
2. Modèles mixtes/hiérarchiques
\[ y = \alpha + \beta x +
\epsilon \]
ou
\[ \mu = \alpha + \beta x
\]
\[ y \sim \text{Normale}(\mu,\sigma^{2}) \]
Modèle linéaire hiérarchique simple
\[ y_{ij} = \alpha + \beta x_{ij}
+ \tau_{i} + \epsilon_{ij}\]
\[ \tau_{i} \sim \text{Normale}(0,\,\sigma_{\tau}^{2}) \]
\[ \epsilon_{ij} \sim \text{Normale}(0,\sigma_{\epsilon}^{2})\ \]
où il y a \(i\) groupes et \(j\) observations par groupe
Autre formulation
\[ \mu_{ij} = \alpha_{i} + \beta
x_{ij}\]
\[ \alpha_{i} \sim \text{Normale}(\alpha,\,\sigma_{\alpha}^{2})\ \]
\[ y_{ij} \sim \text{Normale}(\mu_{ij},\,\sigma_{\epsilon}^{2})\ \]Un modèle est dit hiérarchique lorsque certains de ses paramètres sont définis par une probabilité de distribution.
Simulons un exemple
set.seed(12345)
ng <- 20
nobs <- 10
n <- ng * nobs
id <- gl(ng, nobs)
# valeurs fictives des variables explicatives selon une distribution uniforme
x <- runif(n, 0, 100)
# valeurs des paramètres
a <- 100
b <- 1
# on génère l'effet aléatoire pour chaque groupe
rea <- rep(rnorm(ng, 0, 20), each = nobs)
# prédicteur linéaire
mu <- (a + rea) + b * x
# observations à partir d'une distribution normale avec un écart-type de 5
y <- rnorm(n, mu, 5)
# données
d <- data.frame(y, x)Illustrons le modèle
library(lme4)
library(ggeffects)
m <- lmer(y ~ x + (1|id), data = d)
g <- ggpredict(m, terms = c("x", "id [sample=20]"), type = "random")
plot(g, add.data = TRUE, show_ci = FALSE, colors = rep("black", 20), show_legend = FALSE)
Illustrons les pentes par individus
g <- ggpredict(m, terms = c("x", "id [sample=20]"), type = "random")
plot(g, add.data = TRUE, show_ci = FALSE, colors = sample(colors(), 20))
Retrouvons les paramètres
La variance de l’effet aléatoire était de 20² et la variance résiduelle était de 5².
summary(m)## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (1 | id)
## Data: d
##
## REML criterion at convergence: 1314.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.50905 -0.67270 -0.04039 0.61012 2.58655
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 386.9 19.67
## Residual 25.3 5.03
## Number of obs: 200, groups: id, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 101.35098 4.46512 22.70
## x 0.99676 0.01262 78.97
##
## Correlation of Fixed Effects:
## (Intr)
## x -0.153Récupérons ces effets
fixef(m)## (Intercept) x
## 101.3509752 0.9967622ranef(m)## $id
## (Intercept)
## 1 2.8823702
## 2 -22.7874952
## 3 5.4199340
## 4 -28.1373544
## 5 1.0787346
## 6 -10.6043879
## 7 -6.4207965
## 8 31.3275567
## 9 -10.5628148
## 10 4.4707573
## 11 -24.5843971
## 12 -26.8974001
## 13 24.9063468
## 14 23.4502635
## 15 -2.0989706
## 16 -12.2974161
## 17 -1.4157682
## 18 8.9803887
## 19 44.1496752
## 20 -0.8592262
##
## with conditional variances for "id"apply(ranef(m)$id, 2, sd)## (Intercept)
## 19.60505Coefficients
coef(m)## $id
## (Intercept) x
## 1 104.23335 0.9967622
## 2 78.56348 0.9967622
## 3 106.77091 0.9967622
## 4 73.21362 0.9967622
## 5 102.42971 0.9967622
## 6 90.74659 0.9967622
## 7 94.93018 0.9967622
## 8 132.67853 0.9967622
## 9 90.78816 0.9967622
## 10 105.82173 0.9967622
## 11 76.76658 0.9967622
## 12 74.45358 0.9967622
## 13 126.25732 0.9967622
## 14 124.80124 0.9967622
## 15 99.25200 0.9967622
## 16 89.05356 0.9967622
## 17 99.93521 0.9967622
## 18 110.33136 0.9967622
## 19 145.50065 0.9967622
## 20 100.49175 0.9967622
##
## attr(,"class")
## [1] "coef.mer"Complexifions le modèle
b## [1] 1reb <- rep(rnorm(ng, 0, 0.5), each = nobs)
mu <- (a + rea) + (b + reb) * x
y <- rnorm(n, mu, 5)
d <- data.frame(y, x)Illustrons les pentes
m <- lmer(y ~ x + (x|id), data = d)
g <- ggpredict(m, terms = c("x", "id [sample=20]"), type = "random")
plot(g, add.data = TRUE, show_ci = FALSE, colors = rep("black",20), show_legend = FALSE)
Retrouvons les paramètres
La variance de l’effet aléatoire sur la pente était de (0.5)².
summary(m)## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (x | id)
## Data: d
##
## REML criterion at convergence: 1370.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.19310 -0.58325 0.02874 0.58522 2.42595
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 388.2872 19.7050
## x 0.2105 0.4588 -0.27
## Residual 22.0018 4.6906
## Number of obs: 200, groups: id, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 100.4564 4.4699 22.474
## x 0.8959 0.1033 8.672
##
## Correlation of Fixed Effects:
## (Intr)
## x -0.283
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00784532 (tol = 0.002, component 1)Récupérons ces effets
fixef(m)## (Intercept) x
## 100.4563555 0.8958772ranef(m)## $id
## (Intercept) x
## 1 -4.259743 -0.11207853
## 2 -24.555894 0.09388618
## 3 9.760813 0.08930651
## 4 -23.972706 -0.43249104
## 5 6.379230 0.06853379
## 6 -10.672292 0.92265805
## 7 -4.656865 -0.14388477
## 8 33.097624 -0.22946144
## 9 -7.593899 0.23044679
## 10 -2.310431 0.38486838
## 11 -32.815439 -0.47826203
## 12 -26.583571 0.53612804
## 13 22.768061 0.34741573
## 14 22.280452 -0.58389607
## 15 1.432798 0.75717249
## 16 -8.201353 -0.31288440
## 17 1.063482 0.10990165
## 18 10.101182 0.03586551
## 19 39.745709 -0.88012487
## 20 -1.007156 -0.40309995
##
## with conditional variances for "id"apply(ranef(m)$id, 2, sd)## (Intercept) x
## 19.4359900 0.4557504coef(m)## $id
## (Intercept) x
## 1 96.19661 0.78379868
## 2 75.90046 0.98976339
## 3 110.21717 0.98518372
## 4 76.48365 0.46338617
## 5 106.83559 0.96441100
## 6 89.78406 1.81853525
## 7 95.79949 0.75199243
## 8 133.55398 0.66641577
## 9 92.86246 1.12632400
## 10 98.14592 1.28074559
## 11 67.64092 0.41761518
## 12 73.87278 1.43200525
## 13 123.22442 1.24329294
## 14 122.73681 0.31198114
## 15 101.88915 1.65304970
## 16 92.25500 0.58299281
## 17 101.51984 1.00577886
## 18 110.55754 0.93174272
## 19 140.20206 0.01575234
## 20 99.44920 0.49277726
##
## attr(,"class")
## [1] "coef.mer"Formulations
m <- lmer(y ~ x + (1|id), data = d)
m <- lmer(y ~ x + (x|id), data = d) # équivalent à (1 + x|id)
m <- lmer(y ~ x + (0 + x|id), data = d)
m <- lmer(y ~ x + (1|id) + (0 + x|id), data = d)Corrélations intercepts/pentes
m <- lmer(y ~ x + (1|id) + (0 + x|id), data = d)
Le dernier est équivalent à
m <- lmer(y ~ x + (x||id), data = d)
et signifie qu’on assume que l’intercept et la pente ne
sont pas corrélées ou qu’on ne veut pas estimer cette corrélation. En
pratique, il n’y a pas énormément de raison de faire cela et par défaut
le modèle va estimer cette corrélation.
Effets aléatoires croisés
Supposons que nous avons des individus observés plusieurs fois par année et ce sur plusieurs années.
m <- lmer(y ~ x + (1|id) + (1|year), data = d)Effets aléatoires nichés
Supposons maintenant que les individus sont annuels et qu’ils ne sont
observés que pour des années particulières.
m <- lmer(y ~ x + (1|year/id), data = d)
Équivalent à:
m <- lmer(y ~ x + (1|year) + (1|year:id), data = d)Quand utiliser les effets aléatoires?
Mesures répétées sur certaines unités ou individus
Les observations sont regroupées sous des structures ou ne sont pas indépendantes
Il y a de la réplication à plusieurs niveaux de la hiérarchie
On veut estimer la variance ou généraliser à une population d’unités d’échantillonnage
Formulation matricielle
\[\boldsymbol{y} \sim \text{Normale}(\mathbf{\beta}\mathbf{X}, \mathbf{I}\sigma^2)\]
\[\boldsymbol{y} \sim \text{Normale}(\mathbf{\beta}\mathbf{X}, \mathbf{\Sigma})\]
Effet spatial, temporel, phylogénétique
Autocorrélation spatiale, temporelle, phylogénétique, etc…
Régularisation (shrinkage)
Complete pooling
No pooling
Partial pooling